Mon. Not. R. Astron. Soc. 000, [Tp2] (2013) Printed 16 January 2013 (MN MfcX style file v2.2) 



Chaotic cold accretion onto black holes 



M. Gaspari 1 *, M. Ruszkowski 2 ' 4 , S. Peng Oh 3 

1 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 8574-1 Garching, Germany 
2 Department of Astronomy, University of Michigan, 500 Church Street, Ann Arbor, MI 48109, USA 

3 Department of Physics, University of California, Santa Barbara, CA 93106, USA 

4 The Michigan Center for Theoretical Physics, 3444 Randall Lab, 450 Church Street, Ann Arbor, MI 48109, USA 



16 January 2013 



ABSTRACT 

Bondi theory is often assumed to adequately describe the mode of accretion in as- 
trophysical environments. However, the Bondi flow must be adiabatic, spherically 
symmetric, steady, unperturbed, with constant boundary conditions. Using 3D AMR 
simulations, linking the 50 kpc to the sub-pc scales over the course of 40 Myr, we sys- 
tematically relax the classic assumptions in a typical galaxy hosting a supermassive 
black hole. In the more realistic scenario, where the hot gas is cooling, while heated 
and stirred on large scales, the accretion rate is boosted up to two orders of magnitude 
compared with the Bondi prediction. The cause is the nonlinear growth of thermal in- 
stabilities, leading to the condensation of cold clouds and filaments when t coo \/ts 10. 
The clouds decouple from the hot gas and fall toward the centre. Subsonic turbulence 
of just over 100 km s _1 (M > 0.2) induces the formation of thermal instabilities, even 
in the absence of heating, while in the transonic regime turbulent dissipation inhibits 
their growth (iturbAcooi !)• When heating restores global thermodynamic balance, 
the formation of the multiphase medium is violent, and the mode of accretion is fully 
cold and chaotic. The recurrent collisions, shearing and tidal motions between clouds, 
filaments and the central torus cause a significant reduction of angular momentum, 
boosting accretion. On sub-pc scales the clouds are channelled to the very centre via 
a funnel. A good approximation to the accretion rate is the cooling rate, which can 
be used as subgrid model, physically reproducing the boost factor of 100 required 
by cosmological simulations, while accounting for the frequent fluctuations. Since our 
modelling is fairly general (turbulence/heating due to AGN feedback, galaxy motions, 
mergers, stellar evolution), chaotic cold accretion may be common in many systems, 
such as hot galactic halos, groups, and clusters. In this mode, the black hole can 
quickly react to the state of the entire host galaxy, leading to efficient self-regulated 
feedback and the symbiotic Magorrian relation. Chaotic accretion can generate high- 
velocity clouds and cause strong variations of the AGN luminosity, jet orientation, and 
spin. During phases of overheating, the hot mode becomes the single channel of ac- 
cretion, though strongly suppressed by turbulence. Future high-resolution data could 
determine the current mode of accretion: assuming quiescent feedback, the cold mode 
results in a quasi flat temperature core as opposed to the cuspy hot mode profile. 

Key words: accretion - black hole - hydrodynamics - galaxies: ISM/IGM/ICM - 
instabilities - turbulence - methods: numerical 



1 INTRODUCTION 

Accretion onto compact objects plays central role in the evo- 
lution of astrophysical systems. A number of physical pro- 
cesses, such as radiative cooling, heating and turbulence, 
as well as the presence of magnetic fields and transport 
processes, could all affect the effective accretion rate onto 
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compact objects. The complexity of the real accretion pro- 
cess makes it impossible to predict the accretion rate an- 
alytically. Only under very restrictive assumptions of adi- 
abaticity, spherical symmetry, unperturbed and steady ini- 
tial conditions, in the absence of magnetic fields, self-gravity, 
and feedback, can the nonlinear hydrodynamic equations be 
solved analytically ( |Bondi|[l952[ ). In this idealised case, the 
Bernoulli equation reduces to 
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with gas pressure p, density p, velocity v, adiabatic index 7, 
and mass of the accretor M.. After some algebraic manipu- 
lation, conservation of mass and Eq. lead to the famous 
Bondi accretion rate formula 



Mb = A4tt(GM.) 2 Poo/c 



(2) 



where A = (l/2)^ +1) ^"'- 1) [(5 - ^J/^-Cs-st)/ 3 ^-!) is a 
normalisation factor of order unity, and c s is the gas sound 
speed (the infinity symbol denotes very large radii). 

This elegant and convenient formula for mass accretion 
has been extensively used in the last fifty years in countless 
astrophysical studies. In particular, a large majority of in- 
vestigations, either theoretical/numerical or observational, 
assume that that accretion dynamics follows the Bondi so- 
lution even when the region of influence of the accretor 
(the Bondi radius, re = GM./d? j00 ) is not resolved (e.g. 



Reynolds et al. |1996 Locwcnstcin et al. 2001| |Churazov 


et al. 


2002; Baganoff et al.| 2003| Di Matteo et al.| 2003, 


20051 


Springel et al. 2005 |Allen et al. 2006 Croton et al. 


2006 


Hopkins et al.|2006| Rafferty et al.|2006 |Sijacki et al. 


2007[ 


Hardcastle et al.||2007| |McCarthy et al.||2008| |Booth 


& Schayc 2009; Cattaneo & Tcyssicr 2007, Cattaneo et al. 


2009 


Yang et al. 112012 1. However, since realistic astrophys- 



ical conditions are dramatically different from the classic 
scenario, it is essential to include more realistic physics and 
relax the Bondi assumptions. This is especially important 
in the context of accretion onto supermassive black holes, 
as the SMBH feedback governs the formation and evolution 
of galaxies, groups, and clusters throughout the cosmic time 
( McNamara fc Nulsen|2012"l > . 

In the last decade, investigators have started to quan- 
tify the deviations from the classic Bondi accretion scenario, 



using numerical (e.g. Proga & Begelman 2003 Pen et al 



semi-analytical approaches (e.g. 



Soker 



2011 



[2006{ |Pizzolato fc Soker 



Mathews & Guo 



2003 Krumholz et al.||2005[ |2006| |Igumenshchev||2006[) or 

~ Quataert fc Narayan||2000[ 



[2010[ |Narayan fc Fabian] 
20121. However, the latter models 



are limited by restrictive assumptions on the flow properties 
(e.g. steady, spherically symmetric, isobaric, constant cool- 
ing function, free boundary conditions, simplified potential), 
while the former have often focused on particular physics or 



specific classes of objects, as quasars or supernovae (Barai 
et al.||2012 Hobbs et al.| j,2012). Furthermore, previous sim- 



ulations typically studied either the region within rs for a 
few Bondi times (£b = ?"b/c s ,oo), or the opposite regime in- 
volving the kpc and Gyr scales. Other relevant numerical 
limitations include the choice of the geometry (ID, 2D) or 
the discretisation method (SPlQ. 

In the present work, we intend to provide a coherent 
and unified picture of the departures from the classic and 
idealised Bondi accretion scenario. Using three-dimensional 
adaptive mesh refinement (AMR) simulations, we link the 
small sub-Bondi spatial scales to the large 50 kpc scales, 
and study accretion on timescales of over 40 Myr (> 200 
£b). As reference case, we consider accretion onto a black 



1 SPH codes, often used in the past works, have severe difficulties 
in tracking gas mixing, instabilities and inflow boundaries. 



hole in a common astrophysical environment - an ellipti- 
cal galaxy in a hot gaseous halo. We progressively increase 
the realism of the accretion process, while including physics 
that is simple and generalizable to a wide range of situations. 
First, we start from purely adiabatic evolution in a stratified 
atmosphere (jj3). Second, we relax the assumption of adia- 
baticity and include radiative cooling (Q. Third, we drive 
turbulence to model the effect of weak or strong gas stirring 
(e.g. AGN feedback, galaxy motions, mergers, cosmic flows, 
stellar evolution), studying the dynamics with and without 
cooling (ij5]and ijBJ. Finally, we introduce distributed heat- 
ing to emulate the effects of AGN and stellar feedback (*J7|. 
This allows us to keep the system in global - though not 
local - thermal equilibrium, as observed in most galaxies, 
groups, and clusters. 

After taking into consideration the above realistic pro- 
cesses, the accretion dynamics profoundly differs from the 
classic case. The mode of accretion is cold and chaotic, driven 
by stochastically forming and moving multiphase filaments, 
which condense via nonlinear thermal instability (TI). The 
cold clouds are complex extended structures, as opposed 
to infinitesimally small objects on ballistic orbits. The fre- 
quent inelastic collisions, fragmentations, and shearing mo- 
tions between clouds, filaments, and the central torus, effi- 
ciently reduce the angular momentum of the cold phase. The 
amount of gas that is accreted trough the central region is 
thus boosted, up to 100 times the Bondi rate. On sub-parsec 
scales (down to our highest resolution of few tens gravita- 
tional radii), the cold gas is channeled toward the black hole 
through a funnel, maintaining a similar high level of accre- 
tion (j j7.5[ ). We show that the cooling rate is a good ap- 
proximation to the accretion rate, and we suggest to employ 
it as fiducial subgrid prescription for accretion (and thus 
feedback), especially in cosmological simulations. In |JsJ we 
present an in-depth discussion of the results, along with the 
limitations of our models and the role of additional physics. 



2 PHYSICS & NUMERICS 
2.1 Initial Conditions 

We study accretion physics in a very common astrophysical 
system, i.e. an elliptical galaxy embedded in a hot atmo- 
sphere (intragroup/intracluster medium - IGM/ICM). The 
initial density profile is in hydrostatic equilibrium (neglect- 
ing the black hole) with the temperature profile correspond- 
ing to that observed in the typical massive galaxy/group 
darker line; cf. Gaspari et al.|2011b l. 



la 



GM./ci^ ~ 85 pc 
1 on kpc-scales away 



NGC 5044 (Fig. 

The gravitational potential in the centre is dominated 
by the supermassive black hole of mass M. = 3x 10 9 Mq 
Schwarzschild radius is thus Rs = 2GM,/c 2 ~ 3 x 10~ 4 pc 
while the initial Bondi radius is tb ~ ' ' 1 ' ~ 
(the sound speed is c Si00 — 390 km s 
from the black hole), and the corresponding Bondi time is 
tB = rB/c s ,oo — 210 kyr. Adopting a large black hole mass 
increases the Bondi radius and thus allows us to resolve this 
radius better, though the timestep becomes smaller (due to 
higher central temperature and velocities). 

On kpc scales, the total potential is dominated by the 
deVaucouleurs profile of the cD galaxy with stellar mass 
M* ~ 3.4 x 10 11 M Q and effective radius r c ~ 10 kpc. Within 
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the central ~300 pc, M* is assumed to be constant as com- 
monly observed ( Mathews fc Brighenti|2003 1 , thus inducing 
a flat gas density core. On larger scales (tens of kpc), the 
major contribution to total gravity g comes from the dark 
matter potential approximated via the NFW profile with 
virial mass M v i r — 3.6 x 10 13 Mq (r v i r — 860 kpc) and con- 
centration parameter c ~ 9.5, in the concordance ACDM 
universe (Ho ~ 70 km s _1 Mpc -1 ). The normalisation of 
the gas density is set by the gas fraction / gas ~ 0.1 (at r v i r ) 
yielding central electron number density n Cj o — 0.5 cm -3 . 

The ratio of the cooling time, t coo i = 1.5 nkBT/n e riiA 

({2^, to the free-fall time, t« = (2r/g) 1/2 , starts below 10 
in the range ~ 100 pc - 8 kpc. This is typical of cool-core 
systems which constitute the majority of groups, clusters, 
and massive galaxies. In the case of our galaxy the minimum 
ratio is ~ 5 at 250 pc. The value of 10 corresponds to the 
threshold below which the growth of thermal instabilities 
becomes possible ( Gaspari et al.|2012a McCourt et al.|2012 



Sharma et aL||2012| >. This threshold can also be related to 



2/3 



30 keV 



the observed entropy threshold, K = k^T/n. 
cm 2 , below which substantial star formation is triggered and 



extended multiphase gas is observed (Rafferty et al 
Cavagnolo et al.|2008 1 



2008 



We use a modified version of the AMR FLASH4 code 
( Fryxell et al. |2000[) to integrate the equa tions of hydrody- 
namics (see j^3|and |Gaspari et al.|2011a|b for more details). 
The computational box is ~ 52 kpc on a side with the black 
hole located at the centre. We follow the evolution of the 
gas for 40 Myr. Such long evolutionary times (~200 £b) in 
a domain that large (~600 rs) are necessary in order to un- 
derstand the role of turbulence and cooling on the formation 
of extended cold gas at large radii, and to follow the subse- 
quent accretion onto the central black hole. Large dynamical 
range also removes any undesirable influence of the bound- 
ary conditions on our results. For the sake of completeness, 
we assumed boundaries where the gas is allowed to flow out 
of the computational domain, but is not allowed to enter it. 
The use of three dimensions is crucial to study the stochastic 
dynamics and the instabilities, otherwise leading to spurious 
compression, fragmentation, and/or condensation. 

In order to link the pc/sub-pc scales to the scales of 
tens of kpc over long evolutionary times, our only option is 
to adopt concentric fixed AMR mesh. We typically employ 
14 levels of refinement, reaching up to 21 levels when we 
study accretion on sub-pc scales. This implies an effective 
linear resolution of 65536-8388608 zones, i.e. a dynamical 
range of almost 10 million! This is one of the first studies 
achieving such an extended range (cf . |Levine et al.|2008[ ) . 

Each 'spherical' shell covers over 50 cells in the radial 
direction at each level of refinement. Typical resolution of 
the reference runs is ~0.8 pc. In a few runs we reach the 
sub-pc resolution, down to ~20 Rs ■ At this point additional 
physics certainly becomes important and we intend to ad- 
dress this in future work. As shown by the convergence runs 
(j ]7.5[ ), this resolution is nevertheless sufficient to track the 
key features of chaotic cold accretion. 



2.2 Black Hole Sink 

The main difficulty in properly modelling a sink region lies 
in deciding how much mass passing through the accretion 



sphere (with a radius r a ) should be removed from the com- 
putational domain. If the flow is supersonic, all of the gas 
crossing r a should be removed because the external region 
is causally disconnected from the region within the sphere. 
In an idealised (adiabatic) case, the gas never reaches super- 
sonic speed because for 7 = 5/3 the sonic radius is r s = 0. 
In a subsonic flow it is necessary to remove just enough gas 
to avoid an overpressure bounce of the gas that is being ac- 
creted, while ensuring that the gas removal does not lead to 
central depressions in the pressure distribution. 

One of the best methods we found is a self-regulated 
gues^] based on a modified instantaneous inflow rate in the 
central 1-2 zones: M. — Ar 2 p (c 2 + jj^) 1 ' 2 , where A is a 
normalisation factor of order unity (linked to the numeri- 
cal discretisation, and calibrated via the analytic solution). 
For adiabatic flow, this prescription reproduces very well the 
pure Bondi evolution, even if not all of the accreting gas is 
removed from the domain. However, for non-adiabatic (ra- 
diative) flows, the situation is far more complex (j 4)7 1, and 
finding a self-regulated prescription for the accretion rate is 
difficult because of the presence of extended cold gas in free 
fall, which tends to block very rapidly the accretion sphere. 

The best general and robust sink method results in ap- 
plying a vacuum region, by simply evacuating the gas that 



crosses the accretion sphere. For numerical stability (cf. Ruf- 
fert||1994 1 , the gas density is not exactly set to zero, but is 



given a value of 10 -35 g cm -3 in a radius of 4 - 5 zones. 
Within the sink region, velocities are reset to zero and the 
gravitational potential is softened (though the latter is not 
cruciaQ. Throughout this work, we refer to this sinked gas 
mass per timestep as the 'black hole accretion rate', M,. 

The vacuum region is essential to properly evacuate the 
condensed cold gas, avoiding overpressure bounces or un- 
realistic large accretion rates. On the other hand, the vac- 
uum sphere defines a new (pc-scale) sonic point, which at 
first might seem artificial for black hole accretiorj^] How- 
ever, a more realistic relativistic potential makes this as- 
sumption physical, even in the adiabatic case, for the fol- 
lowing reason. Assuming the relativistic potential 4>pw = 
-GM./(r - Rs) ( |Paczynski_fc Wiita||1980| ), and substitut- 
ing it in the Bernoulli Eq. |l]), yields 



= (5 - 3 7 + 8R S + A 1/2 )/8, 



(3) 



where A = (37-5-8i?s) 2 -64 R s (l-J+Rs); the sonic and 
Schwarzschild radius are normalised to rs (hat symbol). The 
key point is that the physical sonic radius resides now at a 
finite distance from the centre. This means that our simu- 
lated flow, adiabatic or not, always becomes supersonic near 
the parsec region, justifying the complete mass evacuation 
via the vacuum sphere. 

Since the sonic point determines M m , using the vacuum 
accretor or the relativistic potential results in a similar mild 
increase in the accretion rate. The non-dimensional accre- 



2 If the estimate for M, is slightly overshooting compared with 
the real solution, then in the next dt it will slightly undershoot. 

3 For instance, <f> = —GM t / (r 2 + s 2 ) 1 / 2 with a softening param- 
eter s 2 = (0.7r a ) 2 exp[— (r/0.7r a ) 2 ]; we also tested Plummer and 
cubic spline softening and did not find significant differences. 

4 In general, the accretor has a finite hard surface, e.g. white 
dwarfs, neutron stars, planets. 
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tion rate (Eq. [5| for the adiabatic reference case with the 
relativistic potential is in fact given bj|^] 



ln+1 



A 



7 + 1 



n 



i?s 7 + 1 



(4) 



where r t corresponds to the accretor surface, if f a > f s , 
or the sonic point in Eq. j3j, if f a ^ f s . For the adopted 
7 = 5/3, A = [f t /(2(h-Rs)) + 3r t /4] 2 and f s = R s + 
(2i?s/3) 1 ^ 2 , instead of the classic Bondi values A = 1/4 and 
f s = 0. Note that the difference in the accretion rate with 
respect to the classic case is moderate (< 20 per cent, for our 
typical setup) , while the key relativistic correction resides in 
avoiding the singular sonic radius. The more general formula 
given by Eq. Q is useful for detailed comparisons with the 
adiabatic (Bondi-like) reference well as for future 

accretion studies. 



2.3 Radiative Cooling 

The hot plasma in galaxies, groups and clusters emits radi- 
ation mainly in the X-ray band due to Bremsstrahlung at 
T > 10 7 K and line emission at lower temperatures. The 
radiative emissivity of the gas is proportional to the elec- 
tron and ion number densities £ = n riiA(T, Z). As in our 
previous works, the cooling function A is modelled follow- 



ing the work of Sutherland & Dopita ( 1993 1 and assuming 
metallicity of Z ~ 1 Zq (Rasmuss en fc Ponman|2 009). The 
stable cold phase has a temperature floor set at 10 K (which 
can be justified by UV background heating). We adopt for 
simplicity a constanl|^] molecular weight of fi ~ 0.62. 

For numerical integration of the cooling term, we use 



the exact cooling module described in Gaspari et al. ( 2012a) 



Even if not strictly required, we adopt a timestep limiter 
imposed by the cooling time, in order to achieve higher ac- 
curacy. We tested also a few runs with an explicit Runge- 
Kutta cooling solver and found comparable results, albeit 
the explicit solver usually requires smaller timesteps. 



2.4 Turbulence Driving 

We model continuous injection of turbulence via a spec- 
tral forcing scheme that generates statistically stationary 



velocity fields (Eswaran fc Pope 1988 Fisher et al. 20081 



This scheme utilises an Ornstein-Uhlenbeck (OU) random 
process (coloured noise), analogous to Brownian motion in 
a viscous medium. The driven acceleration field is time- 
correlated with a zero-mean and constant root-mean-square 
value. The time-correlation is important for modelling real- 
istic driving forces. In the OU process, the value of the gas 
acceleration at previous timestep a n decays by an exponen- 
tial damping factor / = exp(— dt/rd), where Td is the cor- 
relation time. Simultaneously, a new Gaussian-distributed 
acceleration with variance cr 2 = e* /t c \ is added in the fol- 
lowing way: 



From Eqs 



[^jand^ 



A = [g(f = r\)//(M 



f = M 4 /^+ 1 [l/2 + l/(Af 2 ( 7 -l))], g = 



1)] 2(-t-i) , where 
[l/(f B -R s ) + 

1/(7 — 1)] , and M the Mach number normalised to Coo p 
The sonic radius in Eq. |3| is given by the minimum of g(f). 
° fi starts to significantly increase only below ~ 1.6 X 10 4 K. 



-(7-l)/2 



= fa„ + <T a \A - P G n , 



(5) 



where G n is the Gaussian random variable, e* is the specific 
energy input rate, and a n +i is the updated acceleration. The 
six phases of the stirring modes (three real and three imag- 
inary) are evolved in Fourier space and then converted to 
physical space. In this approach, turbulence can be driven 
by stirring the gas on large scales and letting it cascade to 
smaller scales. This is an efficient approach as the alterna- 
tive would involve executing FFTs for the entire range of 
scales, where the vast majority of modes would have small 
amplitudes. In most runs, we impose a divergence-free con- 
dition on the subsonic turbulent velocity. In §5.4)7.4| we test 
transonic turbulence, including the compressional modes. 

The key physical quantity of interest is the final 3D tur- 
bulent velocity dispersion, a v , which affects the dynamics of 
accretion. The driving of turbulence is intentionally kept 
simple as our goal is not to consider any specific stirring 
source, but rather keep the calculation fairly general. For 
example, the true source of turbulence may be galaxy mo- 
tions, substructure mergers, supernovae or AGN feedback. 
In recent works, Gaspari et al. (2011b 2012b I showed that 



tended periods of time following an AGN outburst. Recent 
observations found very similar values of turbulence in el- 
lipticals, groups, and clusters ( |de Plaa et al . 2012, Sanders 
fc Fabian|20"l2 |. 

In the reference runs the gas is stirred with a v ~ 100- 
180 km s _1 (Mach number M = a v /c s ~ 0.3-0.4). This is 
achieved by adjusting the energy per mod^] e* and corre- 



lation time Td (usually ~ 10 



and 3.15 x 10 11 s, 



respectively). As long as different choices of these parame- 
ters result in the same velocity dispersion, the dynamics of 
the flow remains unaffected. We stir the gas only on large 
scales, L <; 4 kpc, letting turbulence to naturally cascade to 
smaller scales. In the absence of gravity the velocity field is 
approximately Kolmogorov-like, but in a stratified medium, 
as we consider here, the turbulence spectrum is slightly dif- 
ferent. Finally, since the intensity of stirring in our fiducial 
models is low, the turbulent heating, which is proportional 
to tJv/L, is negligibly small (see s|6]). 



2.5 Global Heating 

The final set of more realistic simulations includes heating 
of the gaseous atmosphere. Numerous XMM and Chandra 
observations of the hot gas in galaxies, groups, and clusters 



(e.g. |Vikhlinin et al.|2006| |Diehl fc Statler|2008] |Rasmussen 



& Ponman 2009 ; Su n et al.|2009 1 have clearly demonstrated 
that the majority of systems have a cool core in roughly 
global thermal equilibrium. Despite the presence of radia- 
tive cooling, central temperatures do not drop below one 
third of the virial temperature. Spectroscopic data shows 



that the cooling rates are severely quenched (Peterson & 
Fabian 2006) and the global thermal equilibrium is likely 
maintained at a ~ 10 percent level. The main source of the 



hot gas heating is likely the AGN feedback (Gaspari et al. 
|2011a|b| |2012b| and references within), with other possible 



7 Via simple dimensional analysis jRuszkowski fc Oh 2010j , <r„ 
(jVLe*) 1 '' 3 ; the number of modes is typically N ~ 5000. 
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contributions from mergers, conduction or stellar evolution 
(jBrighenti Sz Mathews 2003). Importantly, these heating 



processes preserve the gentle positive temperature gradient 
seen in the data. We do not model AGN outflows, jets or 
bubbles as in previous work (Gas pari et al.p012c l, but in- 
stead consider a simplified, yet general, heating model which 
encapsulates several phenomena. Specifically, we force the 
heating rate (per unit volume) to be equal to the average 
radiative emissivity in radial shells. In G12b, we showed that 
a prolonged jet feedback heating results in the system set- 
tling roughly to an equilibrium state (i.e. lack of cooling 
catastrophe) characterised by sustained turbulent motions 
mentioned above ({2.4 1. The current heating prescription 



achieves the same goal of keeping the atmosphere in a quasi- 
stable state (c f. |McCourt et al.||2012 i. As demonstrated in 
|Gaspari et al.||2012a| in a jet-heated atmosphere, thermal 
quasi-steady equilibrium is achieved on a timescale longer 
than the duration of the simulations presented here. Since 
now we are not interested in simulating this transition pe- 
riod, we want to ensure that the system settles into global 
thermodynamic equilibrium as quickly as possible. 

The details of how the global equilibrium is achieved do 
not alter the accretion dynamics and our conclusions. Nev- 
ertheless, regarding the shell thickness (where the average 
cooling is computed) , it is preferable to use 4-8 zones in 
order to model the stochastic spatial variation of heating in 
a more realistic way (Fig. 9 in Gaspari et al. 2012b I. The 



very dense cooling gas, T < 10 K, can sometimes produce 
artificial peaks in heating in the inner tiny shells, which can 
be prevented by excluding those outliers from the average of 
emissivity. When stirring is activated, we first let the system 
evolve for a brief period of time in order to seed the pertur- 
bations before imposing global heating. As for turbulence, 
the heating module is fairly general: our findings should thus 
be valid for a range of astrophysical conditions. 



3 ADIABATIC ACCRETION - BONDI 

We start by simulating a purely adiabatic Bondi-like accre- 
tion, i.e. we do not include stirring, heating or cooling. As 
opposed to the classical setup and previous works, we ini- 
tiate the system by using realistic astrophysical conditions, 
i.e. by employing the temperature and density profiles cor- 
responding to a representative galaxy (*J2|. This allows us 
to quantify the differences between the accretion computed 
from the Bondi formula and the numerical simulation that 
are solely due to the non-zero gradients of the thermody- 
namic quantities on large kpc scales. 



3.1 Accretion rate 

The reference accretion rat^] at the initial time is Mb — 
0.09 M Q yr _1 . In the classic adiabatic case (7 = 5/3), p/cf 
is constant (in our stratified atmosphere only up to ~3 re). 
Therefore, the Bondi formula can be conveniently applied 
at ra avoiding the complications due to the gravitational 
potential of the galaxy. 



8 When referring to the reference Bondi formula, Mb, we mean 
the usual Eq. |2J), with the normalisation A provided by Eq. |4j|. 



10 15 20 25 30 
time [Myr] 



15 20 25 30 35 40 
time [Myr] 



Figure la. Adiabatic (Bondi-like) accretion: evolution of the ac- 
cretion rate (1 Myr average). (Top) Accretion rate in physical 
units, Mq yr — 1 . The rate is slightly decreasing due to the pres- 
ence of the galactic gradients. (Bottom) Black hole accretion rate 
normalised to the Bondi formula (j ]2.2| l: the 'boundary conditions' 
are taken at tb (solid) or averaged over 1-2 kpc (dashed). The 
latter, commonly adopted procedure introduces a small bias in 
the accretion rate by a factor of a few (as opposed to ~100 times 
or more that is sometimes assumed in cosmological simulations). 
Note the excellent match between the prediction and the numer- 
ical solution (solid line). 



As shown in Figure 



la 



(bottom), the 3D numerical sim- 
ulation is in excellent agreement with the analytic estimate 
(S 2.2 i, reaching steady state after ~ t^'- the black hole accre- 



tion rate M., through the surface of the sink sphere, is iden- 
tical to Mb (solid line). After circa 10 Myr, the large scale 
gradients start to affect the M. evolution and introduce a 
30 percent decrement in the accretion rate (~ 0.07 Mq yr _1 
at final time; top panel), while the difference with the in- 
stantaneous Bondi rate still remains within few per cent. 
Computing instead the reference Mb on the kpc scale (aver- 
aging over 1-2 kpc) introduces an increase in the normalised 
accretion rate by a factor of 3-4 (dashed line). This small 
bias is again caused by the declining density and slightly in- 
creasing temperature profiles, thus overestimating entropy 
in the Bondi formula. 

Several works adopting subgrid modelling, especially in 
cosmological simulations, boost the Bondi rate by an ad-hoc 
factor a ~100, mainly arguing that the entropy profiles are 
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Figure lb. Adiabatic (Bondi-like) accretion: evolution of the 
mass-weighted electron density (top) and temperature (bottom) 
radial profiles — sampled every 1 Myr, from darker blue to cyan. 
Within ~300 pc from the centre, the profiles are identical to the 
Bondi solution, and they smoothly join with the galactic gradients 
at large radii. The emission-weighted profiles (not shown) are very 
similar. Note the characteristic central cusp in temperature. 



not resolved near re (Booth & Schayc 2009 for a literature 
reviewf^J). However, as indicated by Fig. |la| this bias is typi- 
cally just a factor of a few. On the other hand, we will show 
that the onset of TI and cold accretion ([4|7l can physi- 
cally induce the boosted M, up to two orders of magnitude, 
which seems required to reproduce the observed growth of 
black holes and their host galaxies ( |Di Matteo et al.|2005[ ). 



3.2 Radial profiles &: dynamics 

The radial profiles within 300 pc are identical to the an- 
alytic Bondi solution, with the correct power-law slopes*"] 



at very small radii the density approaches 



-3/2 



and the 



temperature r 



and the slopes of both quantities flatten 



They also note the importance of a multiphase medium, using a 
density-dependent boost factor, which in our cooling simulations 
would result similar toa~ 100. 

10 The limiting slopes are recovered by imposing the free-fall 
velocity vg = (2 GAf./r) 1 / 2 in the continuity equation M = 
iirr' 2 pv, since the gas tends to become transonic at small radii. 
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Figure lc. Adiabatic (Bondi-like) accretion: mid-plane cross sec- 
tion through the mass-weighted temperature at final time (central 
2 kpc 2 ); the normalisation for the quiver velocity field is 500 km 
s _1 . Adiabatic accretion is smooth and steady, even in the pres- 
ence of a galactic potential, and displays a significant temperature 
peak (even when emission- weighted). 



rapidly on larger scales, joining with the galactic/group gra- 
dients (Fig. lb I. The profiles are very smooth and have a low 
scatter, a sign that the code integrates the hydrodynamics 
equations in a consistent and stable way. Over 4 re, the 
profiles are steady and very close to the initial conditions. 

In the case of adiabatic accretion, the emission- weighted 
profiles (not shown) are very similar to the mass-weighted 
ones. This is due to the lack of cold gas. It is important 
to emphasise that the profile of the emission- weighted tem- 
perature T cw within <, re rises substantially. This is due 
to adiabatic compression and is in contrast to the pure un- 
perturbed cooling case that exhibits only a slight decline, or 
the heated case where the profile is roughly flat. With future 
deeper and high-resolution X-ray observations, the mode of 
accretion can thus be unveiled through the measurements of 
the gas temperature profiles. 

The temperature cross section shown in Figure [Tc] il- 
lustrates excellent stability of the simulation: even after 40 
Myr (200 te), there are only tiny oscillations and asymme- 
tries due to the cartesian AMR grids (< 1 per cent). The 
inflow velocity is overall subsonic (in contrast to the cool- 
ing run), approaching M 



1 only near the accretor 02.21. 
As we show below, this spherical (hot) accretion mode, al- 
though very convenient from the analytical point of view, is 
too idealised. Additional physics, such as cooling, stirring, 
and heating, drastically alters the whole evolution. 

We note that several numerical schemes were tested 
(also for the non-adiabatic runs), to check the robustness 
of the results. We tested different flux formulations (split, 
unsplit), Riemann solvers (HLLC, ROE, hybrid), data re- 
construction methods (MUSCL, PPM), characteristic slope 
limiters (minmod, VanLeer, Toro), and key parameters (e.g. 
cfl number, interpolation order). They all give comparable 
results. Overall, we opted for the unsplit formulation with 
PPM plus hybrid solver, which is in principle more accurate. 
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4 ACCRETION WITH COOLING 

Radiative cooling radically changes the evolution of accre- 
tion flows. A fully realistic environment would imply a rough 
balance between cooling and heating Nevertheless, it 

is very instructive to study the accretion flow when adia- 
baticity is violated and cooling prevails. The evolution may 
represent a phase in which the feedback heating is somehow 
delayed (as in the Phoenix cluster; McDonald et al.|2012 l. 



4.1 Accretion & cooling rate 

The accretion rate is very similar to the cooling rate (Figure 
|2a[ ), exponentially increasing until 7 Myr and then saturat- 
ing around 7.5 M Q yr _1 at 40 Myr (~5 icool)- The key result 
is the drastic increase of the accretion rate by > 100 times 
(from 0.07 to 7 M yr _1 ). If we tried to predict M. via the 
Bondi formula, we would be w rong by over 200 times (mid- 



dle panel). As noted in 1 3. 1 even computing poo/c 3 ,oo on 



different scales would change the normalised accretion rate 
just by a factor of a few (see physical M.). This simulation 
demonstrates that the Bondi prescription provides incorrect 
estimates for the accretion rate when cooling is present. Fur- 
ther modifications will be induced by stirring and heating 
(i 5|7 1, but the break of adiabaticity plays the major role. 

It could be tempting to modify the Bondi formula in 
order to retrieve an estimate for the accretion rate. How- 
ever, the latter is suited only for accretion in a single-phase 
symmetric and unperturbed medium. In the current sce- 
nario, the cold gas decouples from the hot flow, and enters 
the quasi free-fall regime (though there is still some fric- 
tion with the hot medium). Therefore, we can infer that 
M. ~ M co \ A /tff (but not Mg aa /iff, especially in the pres- 
ence of TI discussed below). Since the amount of the cold 
gas mass depends on the cooling time, our estimate for the 
accretion rate reduces to M, ~ M e&s /t coo \, which is the cool- 
ing rate. Figure [2a] shows this tight relationship, which will 
still hold in the presence of stirring and thermal instabilities. 
This argument (elaborated in more detail in the following 
Sections) provides a basis for the subgrid modelling based 
on cold feedback, which has been previously shown to be 
extremely efficient in self-regulating the thermodynamical 
evolution of galaxies, groups, and clusters (§81). 



4.2 Radial profiles 

The presence of cooling drastically changes the evolution of 



radial profiles (Figure 2b I. The central flat density profile, 
n c , is erased in less than 1 Myr as in the adiabatic run. 
However, the accumulation of gas is now almost two orders 
of magnitude larger, reaching 100 cm -3 at 10 pc. In addition, 
the density slope is approximately r" 1 ' 5 within two Bondi 
radii (~150 pc; rather than only very near the accretor), and 
it smoothly joins with the galactic gradient with ~ r -1,3 . A 
negative slope of 3/2 is typical of an accreting gas in free-fall 
(j ]3.2[ ). In fact, the cooling gas quickly loses pressure and can 
be thus considered to be in free-fall within few tb- 

The mass-weighted temperature profile (middle panel) 
clearly shows that the whole core has condensed out of the 
hot phase and reached a stable floor at 10 4 K. This is il- 
lustrated by the significant drop in the emission-weighted 
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Figure 2a. Accretion with cooling: evolution of the accretion 
rate (cf. Fig. \la.\ and cooling rate (bottom). The dashed line cor- 
responds to the previous adiabatic run. Note the dramatic boost 
in the accretion rate, by over two orders of magnitude with re- 
spect to the Bondi estimate, and the tight link to the cooling rate. 



profile within 2 re (bottom panel) due to the zero sensitiv- 
ity of X-ray detectors below 0.3 keV. 

One key result concerns compressional heating. Since 
tB is much lower than the central initial t coo \ (~30x), it 
is tempting to study accretion models neglecting cooling 
(e.g 



Narayan & Fabian 20111. However, as the tempera- 



ture profiles show, adiabatic compression can delay cooling 
only for ~ 3 Myr. Even if the ratio of t coo \/te initially in- 
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Figure 2b. Accretion with cooling: mass- and emission-weighted 
radial profiles of density and temperature (cf. Fig. \lb\ . The T ew 
profile has an X-ray cut of 0.3 keV and is computed in larger 
radial bins (to emulate a Chandra observation) . It is evident that 
the massive condensation of cold gas, out of the hot phase, is 
entering the supersonic regime within few re- 



sec region starts to gradually collapscj^] down to 10 4 K. At 
this point, the Bondi solution becomes unrealistic. Adiabatic 
compression can not stop cooling, it can only somewhat de- 
lay the collapse. A quasi adiabatic interstellar medium, even 
within a few parsec from the centre, represents a short tran- 
sient phase in the history of accretion. Interestingly, as we 
explain in the following Sections, the same conclusions hold 
when the spatially distributed heating and turbulence are 
present. These processes also lead to a globally thermally 
stable atmosphere with local thermal instability and the for- 
mation of multiphase gas clouds. 

On the larger galactic scales (r ~ 1-10 kpc), adiabatic 
compression due to the cD potential seems to be able to 
partially inhibit the weak cooling (the T decrement is no 
more than 40 per cent), at least for 40 Myr. However, the 
central temperature dip continues to expand and a massive 
cooling flow develops in a few 100 Myr (with average cool- 
ing rate ~25-30 Mq yr _1 ). Note that the smoothly vary- 
ing emission- weighted T e w profile (with X-ray cut 0.3 keV) 
paradoxically conceals the development of the massive cool- 
ing flow, which is instead exposed by the high cooling (and 
star formation) rates. 



4.3 Dynamics 

The temperature slice in Figure [2c] (left panel) shows com- 
plete condensation of the galactic core into the cold 10 4 K 
phase. The collapse is symmetric and monolithic and no 
thermal instabilities are present. This giant cold core/sphere 
continues to slowly increase in size. This is accompanied by 
a steady increase in cooling rates. The flow becomes super- 
sonic at ~170 pc (40 Myr). The transition to the supersonic 
regime occurs at the boundary of the cold nucleus and the 
hot gas at larger radii. This boundary coincides with the 
region where the gas enters the fast cooling regime < 10 6 K 
(A(T) oc T- 1 ' 2 ). 

The very large sonic radius, compared with that in the 
Bondi problem, is a key characteristic of radiative accretion 
flows (especially in realistic potentials including the black 
hole, galaxy, and dark matter). Even if the initial sonic point 
lied within tb , the mass accretion rate would eventually be- 
come so large pushing the sonic limit r s > ra (cf. |Qua tacrt 
fc Narayan|2000||Li fc Bryan|2012| [Mathews fe Guo|20 12p| 



Nevertheless, cool core systems typically display a minimum 
^cooi/iff - linked to the eventual sonic point - above the 
Bondi radius (e.g. |McCourt et al.||2012[ ). As shown in the 
large-scale density map with the velocity field overlaid (right 
panel), the inflow is significantly enhanced even at the large 
distance of r ~ 20 kpc, which consequently boosts M.. 

We note that the solutions are characterised by good nu- 
merical stability with very low oscillations both in space and 
in time. Achieving this is non-trivial especially in the pres- 
ence of radiative cooling which tends to amplify oscillations. 
The key to prevent such oscillations is to properly model 



creases toward smaller radii (oc r _1//2 , assuming a transonic 
inflow with T oc r _1 ), the dense warm gas near two Bondi 
radii starts to cool substantially creating an exponentially 
colder dip in temperature. After one cooling time (~8 Myr), 
the conditions near one Bondi radius are completely altered 
by radiative cooling, and the temperature in the inner par- 



Note that very low central resolution, as in cosmological simu- 
lations, will create an artificial inversion of the temperature gra- 
dient (VT < 0), in contrast to the findings of Fig. |2b| 
12 From another perspective, the effective 7 is lower than 5/3 in 
the radiative case, which increases the sonic radius. 
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Figure 2c. Accretion with cooling: mid-plane cross sections through T m (central 2 kpc 2 ) and density (52 kpc 2 ) after 40 Myr; the velocity 
field normalisation is 6000 km s~ 1 . The cooling gas becomes supersonic near few r& due to the monolithic and symmetric condensation, 
which dramatically boosts M% . The inflow continues to increase over tens kpc due to the loss of pressure support. 



the sink region (see [2.2 1. The gas must be properly evac- 
uated when central resolution is very high. Otherwise, the 
gas rapidly accumulates and produces severe back-pressure 
on the inflowing gas at larger radii. This modifies the inner 
solution (e.g. v r > 0) and leads to artificial oscillations. This 
is the main reason why it is crucial to adopt the vacuum sink 
region, especially for the case of accretion driven by TI. 

Although adding radiative cooling to hydrodynamics is 
more physical than resorting to the use of the standard adi- 
abatic Bondi solution, the fact that the flow is still perfectly 
symmetric, steady, and unperturbed still implies that this 
simulated accretion does not adequately approximate that 
occurring in a realistic astrophysical system. In the next 
Sections, we increase the realism of the accretion process 
but note in passing that greater realism of these simulations 
does not affect our finding that M, is significantly boosted 
beyond the prediction from the Bondi formula. 



5 ACCRETION WITH COOLING & 
TURBULENCE 

Observed systems are never characterised by perfect spher- 
ical symmetry. In a cooling medium, turbulence can seed 
new structures and nonlinear instabilities in the atmospheres 
provided that the seeded fluctuations are sufficiently large. 
Turbulent motions in the hot medium can be driven by nu- 
merous astrophysical phenomena including AGN feedback, 



galaxy motions, supernovae, and mergers (e.g. Ruszkowski 
|fe Oh|2011||Vazza et al.(20T2} . We keep our model as general 
as possible in order to better understand basic physics of the 
accretion process. Therefore, continuous non-compressive 
turbulence is driven throughout the evolution with a low 
velocity dispersion reaching 150 — 180 km/s when fully de- 
veloped, i.e. M ~ 0.35-0.4. In §5.4| we test stronger turbu- 
lence up to the transonic regime, a less general case, but still 
relevant when the atmosphere is heated by violent events. 



5.1 Dynamics 

For the present case it is more instructive to start by dis- 
cussing cross sections through the temperature and veloc- 
ity distributions (Figure 3a I. The first key result of this 
more realistic model is the formation of dense cold clouds 
(T ~ 10 4 K). This phenomenon is simply explained through 



the growth of the nonlinear thermal instability (cf. Field 
[19651 |Pizzolato fc Soker|[2005l |Sharma et aT] [2012] ) . This 
occurs despite the fact that there is no distributed heating, 
which can promote the growth of even small fluctuations (sJT] 
below) . The moderate-amplitude stirring destroys the initial 
symmetric conditions and induces perturbations in density. 
Overdense regions then start to cool progressively faster. 
The process becomes nonlinear (although not exponentially 
as in the heated case) and cold dense filaments condense 
out of the hot phase after ~10 Myr (Fig. 3a I. As the cold 
gas falls down in the potential well, compression accelerates 
cloud cooling, although in part delayed by adiabatic heating. 

In order to become nonlinear, the amplitude of the den- 
sity fluctuations should be > 10 per cent in the cooling-only 
case. Otherwise perturbations are advected inwards without 



being able to fully develop (e.g. Balbus & Soker 1989 Bin 



ney et al.||2009| |Li fc Bryan||2012| >. The key point is that 
a very common subsonic turbulence (M ~ 0.35) can induce 
sufficient overdensity amplitude to generate TI and extended 
multiphase gas. 

Another crucial result is that the mode of accretion has 
completely changed. Cold accretion is now fully chaotic (due 
to the high sensitivity of the long-term detailed dynamics on 
TI) and stochastic (given the random nature of the fluctua- 
tions^] Cold clouds and filaments form first where i C ooi/te 
has a minimuup 4 ] of ~ 4-5 at 250 pc (see [5.2|. Over time 



the cold blobs and filaments start to appear at progressively 



13 In this work, we alternate the use of 'chaotic' and 'stochas- 
tic', since they are linked (see also ij7]l; note also that specialised 
mathematical work may present different definitions of 'chaotic'. 

14 In the absence of heating, the threshold for multiphase gas 
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Figure 3a. Accretion with cooling &: subsonic turbulence 
(cr v ~ 150 km s — 1 , M ~ 0.35): mid-plane cross section through T m 
(i = 40 Myr). Including moderate turbulence generates perturba- 
tions in the hot medium. These perturbations grow nonlinearly 
via thermal instability and produce cold filamentary structures. 
The dynamics of the gas is chaotic, driven by collisions, shearing 
and tidal motions, which substantially reduce the clouds angular 
momentum and accelerate the sinking rate. 



larger radii, where the ratio t coo \/tg is initially higher and 
where there is less (spherical) geometric compression. The 
accretor region is small compared with the typical cross sec- 
tion of the condensing clouds/filaments, thus they start to 
collide multiple times: in part they are accreted, in part they 
form a transient rotational structure (torus). New clouds fall 
down and collide with the torus, drastically altering its dy- 
namics and shape. The presence of the torus further facili- 
tates accretion, since it increases the cross section for captur- 
ing the clouds. This compensates for the otherwise expected 
decrease in the accretion rate due to the cloud finite angular 
momentum. 

The cold blobs should not be considered tiny bullets 
with insignificant cross section. As we quantify in |5.3[ ac- 
cretion rate is now only slightly reduced compared with the 
unperturbed cooling case, and it is still up to two orders of 
magnitude higher than the Bondi prediction. Overall, accre- 
tion is significantly chaotic rather than ordered with unper- 
turbed ballistic cloud trajectories. 



5.2 Radial profiles 

Unlike in the pure unperturbed cooling run (Q, turbulence- 
induced overdensities lead to the formation of substantial 
amounts of cold gas clumps. This happens within the first 
~5 Myr. Just as in the pure cooling case the temperature 
still drops to ^ 10 5 K, but the stochastic nature of the pro- 
cess is now evident in the mass-weighted profiles (Figure 
3b I. Cold blobs manifest themselves as noticeable density 
fluctuations superimposed on a strong increase in density 
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formation is also proportional to the perturbation amplitude (cf. 
Sharma et al.|2012 i. 



Figure 3b. Accretion with cooling & turbulence (M ~ 0.35): 
mass— and emission-weighted radial profiles of density and tem- 
perature (cf. Fig. |2b| . Note the condensation of warm/cold gas 
out of the hot phase via thermal instability, and the fluctuations 
imparted by turbulence. 



within the central 3 kpc (up to 10 cm ). The average cen- 
tral density is now an order of magnitude higher compared 
with the unperturbed cooling run. 

On the other hand, if the density and temperature is 
weighted with X-ray emission, the picture becomes very dif- 
ferent. The temperature distribution shows less fluctuating 
profiles, and the density is only slightly peaked, up to ~10 



cm (not shown - see Fig. 5b for a similar example). Re 
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markably, T cw (bottom panel) is almost flat, from few pc 
up to 10 kpc, and slightly hotter than the initial condition, 
~10 7 K. We note that high resolution X-ray observations 
could shed light on the mode of accretion by determining 
the temperature gradient in galactic nuclei. Note that the 
mass of the black hole can be no more estimated via the 
temperature peak, as in the adiabatic Bondi-like case. 

The minimum in the £ C ooi/te profile determines where 
the onset of the Tl and cloud formation will begin. In our 
case, this minimum occurs around 250 pc (see Figure [5c| for 
a similar trend). Beyond ~7 kpc, the cooling time becomes 
too long (n c < 10~ 2 cm -2 ) to produce significant cooling 
during the evolution. In the other regime, within the Bondi 
radius 80 pc), the dynamical time is too short for the 
gas to condense while being advected inwards. It is thus not 
surprising that most of the multiphase gas in the core of 
galaxies, groups, and clusters is observed from hundreds pc 
to se veral kpc (|Mathews fc Brighenti |2003| |Rafferty et al 



2008|jMcDonald et al.|2010||2011a||Davis et al.|2012|| Werner 



et al.|2012 I. This circumnuclear region is the nursery of cold 



gas and star formation. In spiral galaxies, where the whole 
system can have t coo \/tg < 10 up to several tens kpc, TI 
clouds may be even more common (e.g. the high-velocity 
clouds in the Milky Way and M31). 

Outside the core (~5 kpc), the profiles are very similar 
to the initial conditions. Only a very small amount of gas 
is flowing out of the box due to turbulence. Therefore, it 
would not be appropriate to study the evolution for several 
tens Myr in smaller domains (<, 10-15 kpc). Furthermore, 
gas outflow from the computational domain would be ex- 
acerbated by too strong turbulence. The large dynamical 
range achieved in the present simulations is thus crucial to 
properly model also large scales. 

Another reason to keep turbulence at low levels is that 
stirring motions heat the gas via turbulent dissipation. The 
specific heating rate due to the dissipation of turbulence can 



be approximated as e tu 



a%jh (oc M , subject to order 



unity uncertainty). Subsonic velocity dispersions of ~150 km 
s _1 [L > 4 kpc) produce negligible heating compared with 
radiative cooling, as shown by the T profiles at large radii 



(Fig. 3b l. On the other hand, strong transonic turbulence 
(j ]5.4[ ) significantly increases the thermal energy of the gas 
and damps thermal instability. 



5.3 Accretion & cooling rate 

The accretion rate shows typical signs of stochastic evolution 
(Figure [3c| . Until 6 Myr the accretion rate is similar to 
the one found in the unperturbed cooling run because the 
cooling is still feeble. At this time, turbulence is already 
well established, with a v ~ 120 km s _1 . After this transient 
period, M, begins to fluctuate and the peaks in accretion 
rate correspond to moments when dense cold clouds cross 
the accretor region. The typical period of fluctuations is ~5 
Myr. 

The accretion rate reaches a peak value of4-5M0 yr~ 1 , 
similar to that found in the unperturbed cooling run at ~30 
Myr. However, the level of accretion is now lower on aver- 
age by a factor of a few. Compared with Bondi theory, the 
average accretion rate is still up to two orders of magnitude 
greater than Mb (middle panel). The relevant result is that 
radiative losses dominate the dynamics of accretion even in 
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Figure 3c. Accretion with cooling & turbulence (M ~ 0.35): 
evolution of the accretion and cooling rate (cf. Fig. |2a| in the 
top panel the rate is averaged over 0.1 Myr). The dashed line 
corresponds to the adiabatic run. Although in the presence of 
turbulent motions (cf. jJSJl, the accretion rate is again boosted by 
up two orders of magnitude with respect to the Bondi rate at a 
given time, showing a tight connection with the cooling rates. 



the presence of moderate turbulence. In the adiabatic case, 
vortical and bulk motions lead instead to a decrement in the 
accretion rate (Sj6j. 

The absence of heating has a strong effect on the cool- 
ing rate (bottom panel), showing values very similar to the 
pure cooling evolution (slightly reduced by turbulent mo- 
tions), with peaks up to 7 Mq yr _1 . Cooling rate also shows 



© 2013 RAS, MNRAS 000,[T}j22] 



12 Gaspari, Ruszkowski & Oh 



stochastic evolution, though there is a steady trend for the 
rate to increase with time. Af. is again linked to the cooling 
rate, with a short delay between gas condensation and its 
infall. 



5.4 Strong turbulence 

We also investigated the effect of stronger stirring on accre- 
tion. This case is perhaps less general but could be applica- 
ble to situations when more vigorous gas stirring is induced 
by strong AGN outbursts, stellar feedback, major mergers. 
Such events could result in transonic/supersonic velocity dis- 
persions and compressible turbulence. Here we describe only 
key similarities and differences between these runs and those 
discussed above. 

First, we considered turbulence with the velocity disper- 
sion of ~300 km s _1 (Af ~ 0.7). This was done by increasing 
the energy per mode e*. As in the previous run (i 5.3 1, we see 
the formation of thermal instabilities and subsequent chaotic 
accretion. Stochastic accretion phase begins after about one 
cooling time since the beginning of the simulation (~7 Myr). 
Accretion rate is again highly fluctuating with slightly lower 
instantaneous peaks reaching 2-3 Mq yr _1 and a similar 5 
Myr fluctuation period. We also tried 100 times longer corre- 
lation time Td, but found that this did not significantly alter 
the evolution. Overall, the results are similar to the previ- 
ous run with cooling and weaker stirring. However, turbulent 
heating starts to become more relevant: away from the black 
hole sphere of influence, on kpc scales, the temperature in- 
creases by ~30-35 percent and the cooling/accretion rates 
are reduced by a similar fraction. 

When the velocity dispersions become transonic the dy- 
namics completely changes. Strong turbulence leads in fact 
to substantial dissipational heating (oc a^), which easily in- 
creases tcooi/iff over 30, damping thermal instability. The 
suppression of the TI growth may be also supported by 
adiabatic processes ( jSanchez-Salcedo et al.|[2002[ ): random 
compressions heat the gas and fluctuations can re-expand 
before having time to cool. However, in the present run with 



<r„ ~ 600-650 km s" 1 (M ~ 1.5) we see significant varia- 
tions in entropy K, meaning that adiabatic processes play a 
secondary role. 

Turbulent dissipation timescale can be approximated 
as tturb ~ eth/eturb ~ Af -2 teddy, where e t h is the spe- 
cific thermal energy and i e dd y — L/a v is the eddy turnover 
time (with a factor of order unity uncertainty). Assuming a 
characteristic scale L ~ 4 kpc and transonic Mach number, 
yields tturb ~ 6-8 Myr (~ tcool), which implies that turbu- 
lent heating becomes very important. In fact, during the 40 
Myr evolution, the temperature raises by a factor of circa 3. 
Cooling is therefore suppressed by two orders of magnitude 
(< 0.01 M yr" 1 ), and so is thermal instability. The evolu- 
tion is dominated only by strong stirring and only hot gas is 
present in the atmosphere. Consequently, Af. shows a pro- 
gressively strong decline, from 0.1 Mq yr _1 down to 10~ 3 
(30 Myr), due to the high temperature and large vorticity 
(see □61) . The minima of Af. display values even one order of 
magnitude smaller than those computed via the Bondi for- 
mula. Strong turbulence tends also to smooth out the central 
density profile and to drive more gas out of the boundaries. 
The decrement in the central density is 1.5- 2x, a sign that 
turbulent diffusion is also quite relevant via the flattening 
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Figure 4a. Adiabatic accretion with subsonic turbulence 
(<T„~150 km s , M ~ 0.35): evolution of the accretion rate 
(cf. Fig. | la) . Note the significant decline in the accretion rate 
compared with the Bondi-like case without turbulence. This de- 
crease is due to the presence of turbulent eddies (vortical and 
bulk motions relative to the accretor). 



of gradients. At these high Mach numbers, dissipation tends 
to dominate over diffusion (with a net increase of entropy), 
albeit the two contributions are not trivial to disentangle. 

In summary, turbulence over 500 - 600 km s" 1 (A/ > 1), 
or more generally when tturb i; tcool, has the exact opposite 
influence on the dynamics than radiative cooling. It dramat- 
ically reduces the black hole accretion rates, while heating 
the gaj^] The only channel of accretion is thus a very feeble 
hot mode. The whole accretion process is overall determined 
by the competition of three timescales: tg, t CO oi and t tUT h- 



6 ADIABATIC ACCRETION WITH 
TURBULENCE 

It is worth to briefly analyse the evolution without radia- 
tive cooling, while still in the presence of moderate levels 
of solenoidal turbulence, a v ~ 150- 180 km s _1 (Af ~ 0.35- 
0.4). The turbulent energy is just ~8 per cent of the thermal 
energy (i5 tU rb — 0.55 Af 2 f? t h). The following evolution could 
represent a phase of slight overheating, e.g. after an AGN 
outburst or major merger, in which thermal instability is 
temporarily suppressed (t coo i/tg > 10) and the hot mode is 
thus the only channel of accretion. In this Section we em- 
phasise only the relevant differences between such a case and 
its radiatively cooling counterpart discussed in Sj5] 

The key result is that the ac cret ion rate declines as op- 
posed to the cooling run (Figure 4a I: Af. drops by a factor 
of a few below 0.07 M yr _1 , which was the minimum value 
reached in the unperturbed adiabatic model. At final time, 
the decrease is a factor of ~2.5. The accretion rate nor- 
malised to the Bondi rate (not shown) is very low: ~2-3 



Mb (the small bias due to the galactic gradients; cf. [3.1 1 



15 Such high turbulence seems to produce too much heating, com- 
pared with observations. Ho wever, a realistically lower level of 
turbulence (200-400 km s~ 1 ; |de Plaa et al.puT^ is an integral 
element of a successful heating model, as it efficiently redistributes 
and isotropises the feedback energy | |Gaspari et ai.|2012b| . 
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Figure 4b. Adiabatic accretion with subsonic turbulence (M ~ 
0.35): mass- weighted profiles of density and temperature (cf. 
Fig. [lb}. The profiles are similar to the Bondi case, though more 
shallow due to turbulent diffusion, and displaying fluctuations. 
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Figure 4c. Adiabatic accretion with subsonic turbulence (M ~ 
0.35): mid-plane temperature cross section (cf. Fig. \lc\ . The 
stochastic nature of turbulence and the driven turbulent eddies 
are evident. 



The M, evolution still shows a stochastic behaviour, though 
with a reduced amplitude. 

When the turbulent motions are significantly subsonic, 
the vorticity driven via the baroclinic instability becomes 
important and inhibits accretion (cf. Krumholz et al 2005, 
2006). This happens because the medium is stratifieof^] On 
the other hand, when turbulence induces a stronger tran- 
sient bulk motion, the reduction of M. is mainly associated 
with the non-zero relative velocity between the accretor and 
the flow motion ( |Bondi fc Hoyle|1944| 



The radial emission- or mass- weighted profiles (Fig. 4b I, 
are very similar to those corresponding to the Bondi-like run 
({ 3.2 I. The only difference is the presence of moderate fluc- 



tuations on the order of few tens of percent, and the ability 
of weak turbulent diffusion to slightly flatten the central 
density profile (which in minor part helps reducing M.). We 
remark that the heating due to the dissipation of turbulent 
motions is negligible for this level of stirring, as indicated by 
the temperature profile at large radii. 

The temperature map (Fig. [4c| reveals the presence of 
turbulent eddies on the scales of a few kpc, i.e. the injection 
scale, which are cascading to smaller scales. The fluctuations 
driven in the current evolution are typically quasi isentropic, 
with significant fluctuations in pressure, as opposed to the 
runs with cooling enabled. As shown in the previous Section, 
increasing a v progressively magnifies the impact of turbulent 
dissipation, hence leading to a mixture of entropy and pres- 
sure fluctuations. Fluctuations with varying entropy prevail 
when turbulence becomes transonic and, of course, when the 
cooling/heating terms are active. 

In summary, adding moderate stirring to the purely adi- 
abatic case drives hot accretion with significantly reduced 
M.. With stronger turbulent motions, the accretion rate is 
reduced by over two orders of magnitude (not shown). In 
less massive galaxies this state should be much more fre- 
quent, since central heating has more severe consequences 
in less bound systems (Gaspari et al. 2012b I, and turbu- 



lent velocities may easily reach M > 1, given the lower gas 
sound speed. SgrA* and NGC 3115 might be two exemplary 
cases | |Wonget al.|2011} . Finally, this run demonstrates that 
the radiative cooling, in combination with turbulence (and 
heating; see below), is the main cause of the dramatic boosts 
of matter channelled to the very centre, compared with the 
Bondi expectation. 



7 ACCRETION WITH HEATING, COOLING & 
TURBULENCE 

Hot gas in galaxies, groups, and clusters is not only cool- 
ing but it also appears to be kept in quasi thermodynam- 
ical equilibrium. The main heating mechanism that pre- 
vents catastrophic cooling is likely AGN feedback, possibly 
helped by mergers, conduction, and stellar evolution. The 
SMBH acts as a thermostat, triggering AGN outflows/jets, 
and subsequent buoyant bubbles, shocks, and turbulent mo- 
tions ( |2.5[ ). The latter phenomena are able to efficiently 
thermalise the kinetic energy, and to isotropically heat the 



16 The process is analogous to the formation of meteorological 
cyclones. 
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medium ( Gaspari et al.|2012a |. The next step in increasing 
the level of realism of our model is thus to include spatially- 
distributed global heating. As our intention is to implement 
physics that is applicable to a range of situations, we im- 
plement a simple prescription for heating that, in agree- 
ment with observations, guarantees global average balance 
of heating and cooling. In this model the gas is globally 
thermally stable, which does not preclude the possibility of 



the development of local TI (£ 2.5 1 . In addition to heating, 
the simulation discussed in this Section includes cooling and 
turbulence at the level identical to that considered in the 
previous Section (a v ~ 150-180 km s _1 , M ~ 0.35), with 
all other key parameters remaining unchanged. Recall that 
for these velocities, the energy injection due to the dissipa- 
tion of turbulent motions is negligible compared with the 
cooling rate. In §7.4| we experiment with stronger stirring. 



7.1 Accretion & cooling rate 

The evolution of the accretion rate reveals the multiphase 



structure of the flow (Figure 5a I. Before a cooling time 7 
Myr), most of the accreted gas is hot (T > 10 6 K): the M. 
evolution is smooth and tends to slowly decline due to the 
driven turbulent motions (cf. |JsJ> , falling below the rate ob- 
served in the unperturbed adiabatic case, M, ~ 0.08 M Q 
yr _1 . From 7 Myr onward, TI becomes exponentially non- 
linear and cold clouds condense out of the hot phase pro- 
gressively boosting the accretion rate. Consequently, the ac- 
creted hot gas remains below 5 per cent of the accreted cold 
phase throughout the rest of the evolution. 

Similarity between the present and the previous 
stirring-plus-cooling run is evident (cf. Fig. |3c| ). However, 
heating introduces important differences. The fluctuations 
in the accretion rate have considerably higher frequency and 
subsequent accretion peaks can appear in just 2 Myr. In fact, 
in a heated atmosphere TI grows exponentially rather than 
linearly, as would be the case in the cooling-plus-turbulence 
case (cf. |McCourt et al.|2012 |. 

Other than driving a more dynamic atmosphere, the 
effect of global heating appears as a significant reduction 
in the average cooling rate, oscillating around 1 Mq yr _1 
(bottom panel), i.e. 15 per cent of the pure cooling rate (in 
the later stage; i 4.1 1. The cooling rate (and M.) has already 



reached a quasi-steady state, while in the pure cooling run it 
continues to rise up to 25-30 Mq yr _1 (after 100 Myr), more 
than an order of magnitude higher. We note that accretion is 
sub-Eddington in all cases considered in this work (MEdd — 
70 Mq yr _1 ), but could be super-Eddington at high redshift 
when the black hole has considerably smaller mass. 

A key result from this experiment is that, even in the 
presence of heating, using the Bondi formula leads to the 
accretion rates underestimated by up to two orders of mag- 
nitude, with typical values around 50 x near final time. This 
is caused by the fact that the dense cold filaments drive ac- 
cretion despite the fact that their volume filling is negligible 
compared with the hot phase. Just as in the cooling case 
with turbulence, the sizes of the clouds and filaments are 
still sufficiently large to facilitate collisions between them. 
This helps to keep the accretion rate at these high levels. 

We remark that our work provides a physical basis for 
the high value of the boost parameter assumed in cosmo- 
logical simulations (a 2> 1), which is required to reproduce 
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Figure 5a. Accretion with heating, cooling & turbulence (M ~ 
0.35): evolution of the accretion and cooling rate (cf. Fig. |2a| in 
the top two panels the rate is averaged over 0.1 Myr). Even in the 
presence of heating, which reduces the average cooling rate, the 
accretion rate is still boosted by TI up to two orders of magni- 
tude, compared with the Bondi prediction. The accreted hot gas 
remains always below 5 per cent of the accreted cold phase. 



the concurrent growth of black holes and their host galaxies 
(e.g. |Sijacki et al.|2007 l. 
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Figure 5b. Accretion with heating, cooling & turbulence (M ~ 0.35): mass- (fop) and emission-weighted (bottom) radial profiles of 
density and temperature (cf. Fig. |2b| . The extended multiphase structure and stochastic nature of the accretion flow is evident within 
10 kpc, but it is largely concealed by weighting with X-ray emission (projection along line of sight will aggravate this effect). The 
emission-weighted temperature profile is remarkably flat: this prediction could be tested with the next generation of X-ray telescopes 
(with resolution < 1 kpc) to discriminate between the hot and cold mode of accretion. 



7.2 Radial profiles 

Since thermal instabilities and filamentary extended cold gas 
grow in the presence of the distributed heating, the mass- 
weighted profiles are dominated again by the dense cold 



phase within few kpc from the centre (Fig. 5b l and show 
sharp increase (decrement) in the density (temperature). As 
for M, , the temporal fluctuations have now higher ampli- 
tude and frequency, a signal that the multiphase gas is now 
more clumpy and filamentary. 

In the linear regime, and for a homogeneous back- 
ground, we would see thermal instabilities forming at all 
distances from the centre as long as t coo \/tfi is sufficiently 
small. However, in the presence of a declining density gradi- 
ent, icooi/tff displays a minimum at ~ 250 pc (Fig. |5c[). In- 
deed, we observe the formation of the first small filamentary 
cold structures near this radius. In the nonlinear regime, 
small-scale perturbations have also the potential to reach 
higher amplitudes (Bur kert fc Lin||2000[ ). After this initial 
phase, we observe the formation of cold clouds with a wide 
range of sizes, from distances of several 10 pc up to ~8 kpc, 
where t C ool/tff remains below 10 (Fig. |5b[ -|5c| the radial av- 
erage masks some of the clouds). This value of the cooling 



to free-fall time ratio is the threshold that approximately 
corresponds to the onset of nonlinear cold gas condensation 
| |Gaspari et al.|2 012a Sh arma et al.|2012 I. Instabilities form 
whenever this ratio falls below the threshold in a certain re- 
gion, hence setting the scale of the local cloud. 

We note that turbulence acts as an effective diffusive 
process, like thermal conduction, preventing most of the TI 
to reach the smallest resolution (and thus helping conver- 
gence; [7.5 1. The usually tiny Field length (19651 - the 



smallest scale of TI collapse set by conduction in the ab- 
sence of stratification - is thus not relevant in this context. 

After 40 Myr, the dropout of the cold gas from the 
hot phase causes the ambient medium to become progres- 
sively more tenuous, leading to t C odl/ts ii 10 at all radii. 
The growth of new cold clouds is inhibited, since in this 
regime buoyancy dominates the dynamics following convec- 
tive stability. At this point, a strong AGN feedback event 
will likely be triggered due to the cold accretion of the pre- 
viously formed clouds, followed by a phase of slight over- 
heating (*j6j. Consequently, the feedback input will rapidly 
decrease, allowing cooling to lower t coo \/tff and to start the 
loop again (Gaspari et al. 2012a for a detailed analysis of 
the self-regulated feedback cycle over several Gyr). 
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Figure 5c. Accretion with heating, cooling &; turbulence (M ~ 
0.35): volume- weighted ratio of the cooling time and free-fall time 
for the hot phase in radial shells (T > 0.1 keV; every 2 Myr). 
The threshold of ~10 marks the point where the condensation of 
extended cold gas becomes possible (100 pc-8 kpc). The dropout 
of cold gas out of the hot phase results in a steady increase of this 
ratio, since f C ool of the more tenuous hot gas becomes longer. 



The presence of heating has the effect to sustain the 
temperature radial profile at large radii. In fact, the T profile 
remains adherent to the initial condition over several kpc 
(~10 7 K), while in the cooling-only unperturbed case it is 
steadily declining. Removing the contribution to the profiles 
from the cold gas via emission weighting results in a roughly 
constant temperature profile (Fig. |5b[ bottom right panel). 
The X-ray emission-weighted density profile demonstrates 
that the temporal fluctuation^ 7 ] are larger compared with 
the fluctuations in the no-heating case, especially within the 
inner 100 pc (bottom left panel). 

We emphasise that current X-ray telescopes - with typ- 
ical resolution > 1 kpc - have severe difficulty in determin- 
ing the mode of accretion via the Tew and n e profiles (cf. 
the magnitude of errors in Wong et al. 20111. Figure 5b 



also shows that extrapolating the observed profiles from kpc 
scales (e.g. Allen et al.|2006 I is risky, since the X-ray profiles 
resemble Bondi hot accretion, while the accretion rates and 
dynamics could be profoundly different (the same is valid 
for the adiabatic case with turbulence in Sj6|. 

7.3 Dynamics: chaotic cold accretion 

The temperature and density maps clearly reveal that accre- 
tion in a realistic astrophysical situation is cold and chaotic. 
This is a crucial result with important implications for the 
evolution of galaxies, groups, and clusters (|8j). 

Figure |5d| shows a representative example of the tem- 
perature (top) and density (bottom) distribution. After few 
cooling times, the evolution reaches equilibrium in the sta- 
tistical sense. Thermal instabilities quickly enter the non- 
linear regime and extended cold filaments drop out of the 
hot phase within <J 8 kpc from the centre. In this region 
icooi/tff ^ 10. 

17 The spatial variations are significantly smeared out by the 
large radial bins, commonly used also in observational analysis. 
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Figure 5d. Accretion with heating, cooling & turbulence (M ~ 
0.35). Both maps show representative snapshots of the temper- 
ature (top) and density (bottom) distribution, with the velocity 
field overlaid. The mode of accretion is cold and chaotic, driven 
by frequent collisions, shearing and tidal motions between the 
cold clouds, extended filaments and the central intermittent torus. 
These strong interactions reduce the angular momentum of the 
cold gas, leading to the boost of M # . 



The denser gas at 10 4 K forms the nucleus of the ex- 
tended cold filaments. The massive nucleus is often sur- 
rounded by layers of cooling gas in which t coo \ /tg ~ 4 - 6 and 
temperatures are approximately a few 10 6 K. This regime is 
highly thermally unstable due to strong line emission. The 
scenario is strikingly similar to what is seen in the multi- 
wavelength observations of cool cores. These observations 
reveal co-spatial presence of filamentary cold/warm gas in 
Ha, FUV, and soft X-ray band ( |McDonald fc Veilleux|2009 



McDonald et al. 2010 2011a I. A small offset can arise in 



regions of strong stirring or when AGN outflows are active. 
Nevertheless, the overall correspondence between X-ray and 
Ha emission suggests that the filaments developed through 
thermal instability, like in our models. 

Figure |5d| shows that multiphase gas forms on several 
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Figure 5e. Accretion with heating, cooling & turbulence (M . 
0.35). Entrop y di stribution of the same map presented in the top 
panel of Fig 



5d 



(A' = fcaT/ne^ 3 ). Most perturbations induced 
by thermal instability are not isentropic, but quasi isobaric. Adi- 
abatic processes have therefore a minor relevance in chaotic cold 
accretion. The medium is strongly biphase, with the low entropy 
phase emerging from the high entropy background. 



scales with various topologies and complex dynamics - a 
trademark of chaotic systems. The cold gas is always dense, 
in rough pressure equilibrium with the hot ambient. The 
combined action of tidal forces and turbulent motions can 
easily stretch smaller spherical clouds (top panel) into longer 
and thinner filaments, extending even 2 kpc in length (bot- 
tom). Magnetic forces could further contribute to the devel- 
opment of filamentary structures in the flow, due to the ac- 
tion of anisotropic thermal conduction ( Sharma et al.|20l0 1 . 

As revealed by Figure [5eJ the temperature fluctuations 
are due to entropy variations (non isentropic), while the 
pressure remains rather smooth (quasi isobaric), in contrast 
to the evolution of the adiabatic case with stirring (*|6|. 
Adiabatic heating or cooling due to turbulent motions has 
therefore a minor relevance in chaotic cold accretion (cf. 
j ]4.2[ ), apart from its effect on seeding TI. The final product 
is a strongly biphase medium, with the low entropy phase 
emerging from the high entropy background. 

Collisions between cold clouds are very common and 
may occur more frequently than once per Myr. The collision 
frequency increases as the clouds get closer to the accretor. 
Consequently, larger clouds are formed due to merging of 
smaller clouds. Tangled filaments increase their cross section 
over time (even up to ~1 kpc 2 ) and experience even more 
encounters and shear with other clouds. The morphology 
of the cold phase continuously changes and, as the clouds 
grow in size, they may become more spherical. However, 
larger clouds can again become filamentary or may even get 
dissociated by violent motions and tidal forces. The evolu- 
tion is fully chaotic. In fact, the accretion rate and orbits of 
the clouds at a given time and position can not be exactly 
predicted (in a way similar to meteorological weather). This 
scenario could also explain the frequent variability of the 
AGN luminosity. 

We emphasise that approximating clouds as compact 



and dense clouds on ballistic trajectories (e.g. in semi- 
analytic treatments or in numerical simulations of isolated 
clouds) would not capture the essential physical processes. 
In idealised atmospheres, a moderate amount of angular 
momentum could prevent accretion. However, cold clouds 
are entities with finite cross section, which can often form 
extended threads. This leads to numerous inelastic colli- 
sions with a significant loss of the gas angular momentum 



(chaotic clouds may even reduce the black hole spin, King 
& Pringlc 2006, and induce fast jet precession). The friction 
between the cold clouds with the hot phase, although minor, 
also helps to reduce the angular momentum of the clouds. 
These processes substantially boost the accretion rate (cf. 
the turbulence-only model in [Rjl. Notice that, during the 
secular evolution, the cold gas will not always display an in- 
flowing kinematics, since it can be uplifted by violent feed- 



back events (outflows, jets, bubbles; Gaspari et al. 2012a I, 
sustained by the turbulent weather. 

An interesting intermittent structure present in the flow 
is the nuclear torus (Fig. 5b and 5d I , i.e. cold gas possessing 



higher angular momentum within the central 100 pc. Under 
idealised conditions this structure will likely develop into 
a perfectly rotating disc (Shakura & Sunyaev 19731. Some 



angular momentum transport process would be needed to 
facilitate accretion in that case. In chaotic cold accretion, 
the momentum of this intermittent disc is instead removed 
by dissipative collisions with new infalling clouds that con- 
tinuously bombard the torus (especially when they have op- 
posite momentum). The presence of a dense obscuring torus 
within 100 pc from the centre, which is supported by an ex- 
tensive literature on AGN physics ( Bianchi et al.|2012 for a 
review), is thus a natural consequence of chaotic cold accre- 
tion. The dense TI clouds could also represent the narrow 
and broad line region of AGN, situated on scales greater or 
less than the torus, respectively. 

Self-gravity is not expected to substantially change the 
rate of accretion or momentum extraction from the torus. 
For self-gravity to be important Toomre Q parameter must 
be < 1, or Afdisc <) (H/R)M m , where Mdisc is the critical disc 
mass and H/R is the disc aspect ratio. Assuming aspect ra- 
tio of ~0.1, Mdisc k, 3x 10 s M . We estimate that our typical 
torus masses are < few 10 7 Mq, hence the main channel of 
angular momentum loss is the collisions of the clouds with 
the torus on chaotic orbits. Similarly, the cold clouds and 
filaments are not massive enough to overcome the external 
potential (BH, stars, and dark matter; cf. Li fc Bryan|2012 |. 
Considering just the BH influence, the tidal disruption ra- 
dius is r t ~ (Af.//9 c i oud ) 1/3 > 1 kpc, for p cloud < 10" 22 g 
cm -3 , much greater than that of any cloud (pcloud declines 
with radius). Only a small fraction of the cloud cores, where 
gas is considerably denser, can be Jeans unstable and self- 
gravitating. We defer the effects of star formation and stellar 
feedback to future high-resolution studies. 



7.4 Strong turbulence 

We now briefly discuss models with cooling and heating, but 
consider higher levels of turbulence by increasing the energy 
per stirring mode. 

The first run reaches steady turbulent velocities of ~300 
km s _1 (M ~ 0.7), i.e. about two times larger velocity dis- 
persion than in the previous Section. The overall evolution 
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results to be similar. The accretion rate initially drops to 
lower levels (0.02 Mq yr -1 ) and for a longer time. However, 
after 12 Myr, M. is again boosted by the sinking of cold 
clouds that condense out of the hot phase via TI. The ac- 
cretion rate settles to a quasi-stable stochastic regime, but 
exhibits slightly higher peaks up to 1-2 Mq yr -1 , showing 
shorter period (less than 2 Myr). Similar to the weaker stir- 
ring case, the normalised accretion rate stays in the range 
of 50- 100 Mb, and the cooling rate oscillates around 1 Mq 
yr -1 . The maps show the same stochastic dynamics as found 
in the previous Section, with several collisions, merging, and 
shearing motions. In contrast to the reference case, the ra- 
dial profiles now reveal that turbulent heating starts to be 
effective, increasing the temperature by ~ 40 per cent. 

We also tested the transonic stirring regime, driving 
turbulence up to ~550 km s" 1 (M ~ 1.3), this time only 
for the initial 5 Myr (continuous turbulence suppresses TI; 
[ 5.4 1. This level of turbulence mimics the dynamics of a very 
strong AGN feedback event (~10 45 erg s _1 ). Half the evo- 
lution (15-20 Myr) is dominated by turbulent dissipation, 



showing a drastic decline of the accretion rate (as in [5.4 1 
and cooling rates reduced well below 1 Mq yr . The den- 
sity and temperature fluctuations in the mass- and emission- 
weighted profiles are substantial, showing a T increase by a 
factor of ~2. However, since the stirring is not continuous, 
after this hot accretion period, we observe again the for- 
mation of extended multiphase gas, with M. boosted up to 
circa 1 Mq yr -1 or 100 times Mb- 

We conclude that continuous stirring, in a heated atmo- 
sphere, and for the common Mach numbers M ~ 0.2-0.7, 
leads to the formation of TI and a strong boost in the accre- 
tion rate, while not violating the observational constraints 
on the large-scale profiles of the galaxy/group. On the other 
hand, for transonic stirring, turbulent dissipation becomes 
significant and can stifle TI by increasing the i CO oi/iff unless 
the driving vanishes after few Myr (as in a self-regulated 
feedback loop). The secular accretion history is thus driven 
by the competition of ts versus t coo i versus it ur t>, alternat- 
ing between the cold and hot accretion mode (with the latter 
significantly less efficient). 



7.5 Convergence & the sub-parsec region 

As a last step, we test convergence of the fiducial heated 



simulation with moderate stirring (S 7.1 1 by performing three 



dedicated runs. We focus on the convergence of the accretion 
rate rather than the detailed cloud morphology. 

In the first test, we start the evolution from the same 
initial conditions as before, but we increase the central res- 
olution by 4 times to ~0.2 pc. The system is integrated for 
half the usual evolution time, i.e. until ~ 20 Myr. As shown 
in Figure [fl] the higher resolution run (magenta) is consis- 
tent with the previous reference accretion model (black) 
The cooling rate is also convergent throughout the evolu- 
tion, and radial profiles and maps of density/temperature 
show identical values as in Figure |5b| 

In the second test, the radial extent of the fixed AMR 




Figure 6. Convergence test with 4x higher (magenta) and 2x 
lower resolution (red) than the fiducial model with heating, cool- 
ing and moderate turbulence (black - |7.1} . The accretion rate 
shows that the fiducial run with resolution of ~0.8 pc is well in 
the convergence limit (for the magenta line in the statistical sense, 
since the random noise is not identical). Increasing or reducing 
the resolution leads to a very similar evolution and dynamics. 



shells is halved, thus reducing the resolution by a factor 
of two at each radius. Since the maximum resolution is the 



same (as in [ 7.1 1, we are able to drive stirring with the same 



random numbers in real/phase space. The key result is that 
convergence is clearly achieved (Fig. [fi] red line) , since the 
accretion rate matches that for the higher resolution case 
(black line). Occasionally, the accretion rate in the lower 
resolution case is higher by a few percent due to the for- 
mation of cold blobs over few grid zones. However, most of 
the cold blobs are typically described by a thousand or more 
cells (in volume). Turbulent diffusion acts also as an effec- 
tive diffusion, preventing the cold clumps to collapse to the 
minimum resolution (the details of the cold phase will be 
studied in a following work). We note that while the clouds 
emerge in the entire radial range where icool/tff 10, the 
majority of the them form near the inner better resolved 
radius where this condition holds (Fig. [5c| . The evolution 
of the cooling rate is also convergent, smoothly increasing 
to 1 M Q yr -1 at 20 Myr. Finally, over 40 Myr, profiles and 
maps show identical behaviour described in ^7| We conclude 
that the pc resolution of the fiducial models ensures that the 
results are robust. 

The last test is the most expensive, pushing the 
AMR capability of the code and reaching the resolution 
of ~ 20 Rs ~ 6 x 10~ 3 pc (the dynamical range is almost 
10 million). The time evolution needs to be drastically 
decreased down to 1 Myr. We restarted the fiducial run 
from the final time (40 Myr), keeping the same conditions 
in the 52 kpc box but zooming in on the central region. The 
previous sink region is filled with hot diffuse gas at ~10 7 K. 
After the initial transient phase due to the restarting^] the 



18 Convergence is meant here in the statistical sense as the set 
of random numbers for the Fourier stirring modes is different for 
different resolutions (cr v remains the same). 



19 The system experiences a brief phase of purely hot accretion 
(few 10 kyr). Refilling the older sink region with dense cold gas, 
induces instead an initially higher M m , but the later evolution is 
then identical to the previous case with hot gas refilling. 
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Figure 7. Cross section of the gas density distribution from the 
deep zoom-in test with resolution of ~20 Rg. The accretion rate is 
still maintained at high levels, since the cold clouds are stretched 
and channelled via a funnel toward the BH by its strong gravity. 



multiphase gas steadily accretes reaching an average value 
slightly less than that of the fiducial run in 1 Myr, ~0.4 Mq 
yr _1 . The accretion rate remains high even on sub-pc scales, 



since the strong gravitational attraction (\ 2.2 1 channels 



the cold gas towards the black hole via a narrow funnel 
(Figure [7|. In a few dynamical times, most of the cold 
clouds experiencing chaotic collisions can be potentially 
captured and channeled toward the sub-pc region. Below 
a few tens of Rs, the accretion process is expected to be 
more complex. We will perform dedicated simulations to 
study this region including additional physical processes 
discussed in !|8] 



8 SUMMARY & DISCUSSION 

The simplicity of classic accretion theory makes it extremely 
appealing to apply the Bondi formula (Eq. [2| to the inter- 
pretation of observations, build analytical models, and pa- 
rameterise the black hole accretion rate in large-scale sim- 
ulations. However, the Bondi model requires a number of 
strong assumptions that can be grossly violated in a realis- 
tic astrophysical environment. The Bondi flow must be adi- 
abatic, steady, unperturbed, with no vorticity or heating, 
set by constant boundary conditions at large radii, among 
the most notable. In an attempt to increase the realism of 
the accretion, we presented a systematic investigation of the 
effects of radiative cooling, global heating, and turbulence 
on the black hole mass accretion rates in extreme dynami- 
cal range simulations, and compared the simulated evolution 
with the prediction of the classic Bondi model. Our 3D sim- 
ulations also allow us to isolate the effect of non-vanishing 
gradients in density and temperature on the accretion rate. 
Our model is fairly general and applicable to a range of sys- 
tems, such as hot halos in galaxies, groups, and clusters. 
Below we present the summary of our key findings, discuss 
crucial implications and limitations of the model, as well as 
future directions. 



1. Adiabatic & unperturbed accretion. 

First, we studied adiabatic accretion in the presence of a 
realistic declining density and slightly increasing tempera- 
ture profiles. That is, we relaxed the assumption of the flat 
distribution of thermodynamic variables assumed in the 
Bondi model. When evaluated near rs, the Bondi formula 
provides an excellent estimate for the actual accretion rate 
(within a few percent), even if the accretion rate slowly 
declines over time. Computing instead the Bondi rate on 
the kpc scales (as often done in literature), leads to the 
underestimate in M. by a factor of 2-4. This is due to 
the presence of non-vanishing density and temperature 
gradients; this bias can not reach a factor of > 100. We 
emphasise that the relativistic potential induces a finite 
sonic point (near the pc scale), even in the adiabatic flow. 

2. Cooling & unperturbed accretion. 

The standard Bondi prescription may result in unrealistic 
accretion rates. The main reason for this is the violation of 
adiabaticity: astrophysical hot plasmas are radiating, losing 
thermal energy and pressure support via condensation. In 
a hot radiative atmosphere, the accretion rate is boosted 
by two orders of magnitude, compared with the pure Bondi 
values. After one cooling time, the conditions near the 
Bondi radius are completely altered by cooling, i.e. from 
that point on, accretion is dominated by the thermally 
unstable cold gas rather than the spatially uniform hot 
phase. Within a few 100 pc, the mass-weighted temperature 
profiles drop to 10 4 K and the densities reach over 100 
cm -3 . On the kpc scale, the T profile is only slightly 
decreasing, partially sustained by compressional heating. 
However, after 40 Myr, a strong cooling flow progressively 
alters also these regions (if no feedback is present). 

3. Cooling & turbulent accretion. 

Adding stochastic stirring motions characterised by 
cr„ ~ 100-300 km s" 1 (M ~ 0.25 - 0.7), further increases 
the realism of the simulation. These motions can be 
produced by AGN outflows, bubbles, supernovae, mergers 
or cosmic flows. The key result is the nonlinear growth 
of thermal instabilities and the formation of multiphase 
gas even for moderate turbulence. Accretion is no longer 
symmetric and is instead completely chaotic. Cooling rates 
are still high (> 8 Mq yr" 1 ), and the profiles differ from 
the previous case only for the continuous fluctuations. At 
transonic velocities M > 1, turbulent heating becomes 
significant and starts inhibiting TI and cooling: the key 
threshold is tturb/tcool 1, where ttm-b ~ M~ 2 L/o~ v . 

4. Adiabatic & turbulent accretion. 

Without cooling, the effect of the reference subsonic turbu- 
lence (M ~ 0.35) is instead to reduce the accretion rate by 
a factor of ~3. Accretion is reduced by orders of magnitude 
when M > 1. The profiles resemble the Bondi hot mode 
solution with flatter density and temperature profiles, due 
to turbulent diffusion. In the subsonic regime and in a 
stratified atmosphere, the generation of vorticity via the 
baroclinic instability plays a key role in the reduction of 
M m . When turbulence leads to a transient strong bulk 
motion, the non-zero relative velocity between the accretor 
and the flow further reduces the accretion rate. 
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5. Heating, cooling & turbulent accretion. 
The last and most realistic scenario includes spatially- 
distributed heating. This is motivated by the fact that atmo- 
spheres appear to remain in global thermal equilibrium due 
to self-regulated AGN feedback, possibly helped by conduc- 
tion, mergers and stellar evolution. The heating is broadly 
applicable to a range of situations. The heating offsets cool- 
ing globally, but not locally, leading to the growth of non- 
linear thermal instability. In the heated atmosphere, the 
total cold mass and cooling rate are reduced by an order 
of magnitude. On the other hand, (quasi isobaric) fluctu- 
ations are amplified and TI now grows exponentially. The 
result is a very frequent condensation of cold clouds and 
filaments, in the radial range of ~100 pc-10 kpc, where 
icooi/te Si 10. The key point is that in a common astro- 
physical scenario, accretion is mainly cold and chaotic. The 
dynamics is driven by stochastic dissipative collisions and 
shearing motions which cause a significant reduction of an- 
gular momentum in the cold phase and, hence, a significant 
boost of M,. The accretion rate is fluctuating with peaks 
over 1 M Q yr _1 . Utilising the Bondi formula with the den- 
sity and temperature evaluated on kpc scales would under- 
estimate the true accretion rate by up to a factor of 50 - 100. 

It is remarkable that subsonic turbulence (100-300 km 
s _1 ) can induce thermal instabilities, multiphase gas, and 
chaotic cold accretion. In fact, such a level of velocity fluc- 
tuations seems to be very common in hot gas halos (e.g. 
Schuecker et al.||2004| |Churazov et~aT]|20"T2l |de Plaa et aL] 
2012| ISanders & Fabian 2012 and even higher in the inter- 



stellar medium: Vazquez-Semadeni 2012 1 and may be fur- 



ther confirmed by the upcoming Astro-H mission. 

The cold phase shows a plethora of morphologies on dif- 
ferent scales: from clouds to thin long filaments, stretched 
by the turbulent and tidal motions. A central torus or rotat- 
ing structure is often present, even if its angular momentum 
is stochastically altered by the infalling cold gas. The cold 
phase, although almost in free-fall, can not be considered 
collisionless. A significant cross section of clouds is a crucial 
factor facilitating frequent interactions (especially within 5 - 
10 tb), and in boosting M, beyond the Bondi prediction. 
Collisions of the infalling clouds with the torus also help to 
remove angular momentum from this structure and drive the 
accretion of cold gas. Chaotic clouds may even substantially 
reduce the spin of the black hole (King & Pringlc 2006), in- 
duce fast jet precession, and explain the common variability 
of the AGN luminosity. The ubiquitous presence of cold and 
molecular gas observed in the cores of galaxies, groups, and 
clusters (e.g.|Edge|200"l]|Mathews fc Brighenti|2003||Salome| 
fc Combes|2004||Rafferty et al.|2008||McDonald et al.|2010| 



2011b||Davis et al.||2012| |Russell et al.||2012| |Werner et al.| 
2012 ) seems to corroborate the cold accretion scenario. The 



dense TI clouds could also represent the narrow and broad 
line region of AGN, situated on scales greater or less than 
the torus, respectively. 

The stochastic nature of the accretion process is also re- 
vealed via fluctuations in the radial mass-weighted profiles: 
below 1 kpc the density can exceed 10 3 cm -3 (X-ray bright- 
ness increase would be visible only within ra). The X-ray 
temperature is instead remarkably flat, even down to tens of 
pc. The next generation of deeper observations could deter- 
mine the current mode of accretion through the observations 
of the temperature profiles. The adiabatic Bondi-like model 



predicts in fact a rising X-ray temperature below 0.5 ra, al- 
though this could be also achieved when the central feedback 
becomes very active. Note that extrapolating the (emission- 
weighted) profiles from the kpc scale is instead misleading, 
since cold and hot accretion display the same behaviour. 

We performed convergence tests and obtained consis- 
tent results. Lowering or increasing the resolution, does not 
change the mass accretion rates: M, stays on high levels, 
with similar cooling rates and chaotic dynamics. On sub-pc 
scales the cold phase is stretched and channeled toward the 
accretor via a funnel, while the accretion rate remains on a 
comparable level to that on larger scales. 



8.1 Improvements & additional physics 

In future, we intend to study in more detail the extended 
multiphase gas. This phase not only provides the fuel for 
accretion, but it is also associated with star formation and 
is co-spatial with dust and molecular clouds visible in the 
optical, infrared, and radio bands. Limited resolution may 
favour the condensation of the less resolved clouds and fil- 
aments, which could be also achieved by magnetic draping. 
On the other hand, numerical diffusion tends to prevent the 
initial development of instabilities on the smallest scales, an 
effect that could be mimicked by thermal conduction. All 
these effects can affect the distribution of the detailed emis- 
sion measure as a function of T, which we will study next. 

The magnetic field could play a role in shaping the cold 
phase, while its average low strength of few fiG might be in- 
sufficient to modify the global accretion dynamics. The cold 
filaments could be though characterised by enhanced mag- 
netic fields with vectors aligned along the filaments. Mag- 
netic fields will also play a role in determining the degree 
of mixing between the hot and cold phase, with anisotropic 
thermal conduction likely equalising T along the filaments 
| |Sharma et al.pOlO l. On the other hand, turbulence could 
act as the effective (isotropic) diffusive process, making con- 
duction and the typically tiny Field length secondary. 

Radiation processes may also influence the evolution. 
In particular, Compton heating could have the major im- 
pact, especially near a few Rs, where the cold clouds reach 
extreme column densities and opacities. The coupling with 
the hot medium seems instead inefficient due to its high 
ionisation and low density. A (rare) quasar-like Eddington 
regime is probably necessary to produce sufficient radiation 
pressure to alter the global dynamics. 

In this study, we did not include the strong initial phase 
of BH feedback, focusing more on the fuelling stage. Bipolar 
kinetic outflows or jets can entrain and uplift a fraction of 
the infalling clouds along the direction perpendicular to the 
torus ( Gaspari et al.|[2012a l, and induce higher turbulence 
and new cold phase interactions. Some analytic theories (e.g. 
Blandford & Begelman 1999[) also predict that the gas reach- 



ing the very centre may become unbound. However, the fate 
of the infalling gas depends on the mode of accretion. The 
presence of cold gas condensation near ra means that the 
outer boundary conditions for the sub-pc region would be 
different than normally assumed, and analytic accretion the- 
ories (as ADAF or thin disc) might lead to different results. 
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8.2 Implications & applications 

Considering the chaotic nature of cold accretion, it is not 
trivial to predict the exact accretion rate. However, we find 
that the cooling rates are tightly linked to the accretion 
rates, M, « M coo i, both in the heated and pure cooling 
cases. Due to low resolution, cosmological simulations of- 
ten assume substantially boosted Bondi rates by a factor of 



jacki et al. 



100 (e.g. Springel et al. 2005 Pi Matteo et al. 2005 Si- 



2007 Booth & Schaye 2009), or high feedback 



efficiencies, in order to explain the growth of black holes and 
their host galaxy. Interestingly, we find from first principles, 
that the boost factor is about two orders of magnitude just 
as required by the cosmological models. We suggest thus 
to apply the subgrid prescription M, ai M coo \, which can 
physically reproduce both the increased magnitude and the 
realistic fluctuating evolution of the accretion rate, avoiding 
ad-hoc factors dependent on numeric^] Overall, it is the 
cold mode that drives the BH growth/feedback, explaining 
why large-scale simulations are forced to strongly modify 
Bondi 'hot' accretion in order to match observations. 

The thermal instability scenario provides a direct link 
between the large scales of the hot galactic halo and the rel- 
atively tiny SMBH. The TI criterion determines how much 
material can decouple from the hot halo and rain down onto 
the black hole leading to its growth. The presence of the 
instability and the rate of its development are on the other 
hand controlled by the feedback from the BH. In the cold 
mode, the response time of the black hole to the increased 
rate of accretion (a few iff) is much faster than in the case 
of the hot mode. Therefore, a very efficient feedback can 
be established, leading to the proper thermodynamical self- 
regulation of the galaxy, group or cluster, with a duty cycle 
on the order of the cooling time. 

In this symbiotic mechanism, the Magorrian relation 
(1998) appears naturally, since the accretion rate is po- 



tentially linked to the gas mass of the entire galaxy/bulge 
(tcooi/te 10), while in Bondi hot accretion the supply is 
limited by the region < re- Chaotic cold accretion can thus 
lead to the linear scaling M, ~ M gaSi bui g o ~ /-M^buige- The 
normalisation / depends on the gas-to-stellar mass ratio, 
which is set by the interplay of condensation and feedbaclj^] 
at different epochs; a typical range is / ~ 10~ 3 -10~ 2 , inline 
with the observed scaling relation (McConn ell Sz Ma|2012[ ). 
If stars were collisional systems, the black hole would easily 
swallow them, like the cold clouds, imposing a higher / ~ 1 
(provided adequate external resupply). More massive black 
holes thus reside in more massive galaxies, due to the larger 
available fuel up to few 10 kpc (M. w M CO ol oc M gaSiCoro ), 
rather than due to the larger radius of influence (as in Bondi 
theory where M, oc r B oc M 2 ). This symbiotic link between 
the tiny black hole and the extended environment should find 
its strongest manifestation in the cores of clus ters, where 
black holes can grow up to several 10 10 Mq (McConnell 
et al.|20TT Fabian et al.|20f3 showing also high / values) 



Cold accretion should play a crucial role in high-z galax- 



Diffcrcnt boost parameters lead to profoundly different star 



formation history and black hole growth (Booth &: Schaye|2009||. 
21 Momentum- or energy-driven feedback also determines the 
conversion of M, — M„ to the scaling with the stellar vel ocity 
dispersion, M. oc tr* or oc a%, respectively (cf. Fabian|2012 I. 



ies, since the intergalactic medium appears strongly multi- 
phase, as evidenced by the Lyman-a forest, and the rich gas 
supply has not yet massively converted into (collisionless) 
stars. The massive TI condensation should quickly drive 
accretion to or beyond the Eddington limit, activating the 
quasar regime. At lower redshift, when halos are already well 
formed, the feedback driven by thermal instability should 
become less violent and more recurrent, maintaining the core 
of the system in quasi-stable equilibrium for several Gyr. 

In closing, we remark that the hot mode is not ab- 
sent from the secular evolution. After strong feedback or 
dynamic events, the system typically experiences a phase 
of overheating and substantial turbulence, which is more 
prolonged in less massive/more diffuse halos (isolated ellip- 
ticals, spiral galaxies, etc.), due to the lower binding energy 
and gas sound speed. During this stage, thermal instability 
is suppressed. The only channel of accretion is hot gas, with 
accretion rates severely reduced by turbulent motions and 
heating (^4][6|. NGC 3115 and Sgr A* might be two exem- 
plary cases. Overall, the secular accretion history results to 
be substantially polarised, showing deviations of orders of 
magnitude either above or below the Bondi rate, driven by 
the competition between three key timescales: tff versus t coo \ 
versus iturb (stratification vs. net cooling vs. turbulence). 
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